#for flowchart
library(DiagrammeR)
#to export flowchart as .png
#webshot::install_phantomjs()
library(DiagrammeRsvg)
library(rsvg)
library(grid)
library(xtable)
library(dplyr)
library(ggplot2)
library(gridExtra)
#nest/unnest
library(tidyr)
#map function (kind of like a for loop)
library(purrr)
#tidy model summary
library(broom)
library(readxl)
#for tableby
library(arsenal)

Define geom_split_violin()

Based on https://debruine.github.io/posts/plot-comparison/

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, draw_group = function(self,
    data, ..., draw_quantiles = NULL) {
    data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth *
        (xmax - x))
    grp <- data[1, "group"]
    newdata <- plyr::arrange(transform(data, x = if (grp%%2 == 1)
        xminv else xmaxv), if (grp%%2 == 1)
        y else -y)
    newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1,
        ])
    newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
    if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
        stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <= 1))
        quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
        aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")),
            drop = FALSE]
        aesthetics$alpha <- rep(1, nrow(quantiles))
        both <- cbind(quantiles, aesthetics)
        quantile_grob <- GeomPath$draw_panel(both, ...)
        ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata,
            ...), quantile_grob))
    } else {
        ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
    }
})

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity",
    ..., draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, show.legend = NA,
    inherit.aes = TRUE) {
    layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, position = position,
        show.legend = show.legend, inherit.aes = inherit.aes, params = list(trim = trim,
            scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

Load Data

load('Data/DataWithPropensities_seed1.RData')

#convert PrimaryDiagnosis to factor
dat3$PrimaryDiagnosis <- factor(dat3$PrimaryDiagnosis, levels = c("Autism", "None"))
levels(dat3$PrimaryDiagnosis)
## [1] "Autism" "None"
dat3$PrimaryDiagnosis <- relevel(dat3$PrimaryDiagnosis, "None")
levels(dat3$PrimaryDiagnosis) = c("TD", "ASD")

Reshape data to combine motion QC levels

# create dummy factor to include all subjects
dat3$noExclusion <- ifelse(dat3$ID > 0, "Pass", "Pass")
dat3$noExclusion <- factor(dat3$noExclusion)

# convert KKI_criteria to factor level 'Pass'
dat3$KKI_criteria <- factor(dat3$KKI_criteria, levels = c("Pass", "Fail"))

# convert Ciric_length to factor with reference level 'Pass' to match
# KKI_criteria
dat3$Ciric_length <- factor(dat3$Ciric_length, levels = c("Pass", "Fail"))

# combine Ciric_length, KKI, and noExclusion exclusion into one variable
allVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan", "Ciric_length", "KKI_criteria",
    "noExclusion", "PANESS.TotalOverflowNotAccountingForAge", "SRS.Score", "WISC.GAI",
    "DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw", "ADOS.Comparable.Total",
    "CurrentlyOnStimulants", "HeadCoil", "Sex", "Diagnosis4Levels", "ADHD_Secondary",
    "SES.Family", "Race2", "handedness", "CompletePredictorCases", "YearOfScan",
    "MeanFramewiseDisplacement.KKI")

idVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan", "PANESS.TotalOverflowNotAccountingForAge",
    "SRS.Score", "WISC.GAI", "DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw",
    "ADOS.Comparable.Total", "CurrentlyOnStimulants", "HeadCoil", "Sex", "Diagnosis4Levels",
    "ADHD_Secondary", "SES.Family", "Race2", "handedness", "CompletePredictorCases",
    "YearOfScan", "MeanFramewiseDisplacement.KKI")

qcMelt <- reshape2::melt(dat3[, allVariables], id.vars = names(dat3)[which(names(dat3) %in%
    idVariables)], variable.name = "Motion.Exclusion.Level", value.name = "Included")
## Warning: attributes are not identical across measure variables; they will be
## dropped
# rename exclusion levels NOTE: need None to be highest level for
# geom_split_violin
levels(qcMelt$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")

# convert Included to factor with pass as reference
qcMelt$Included <- factor(qcMelt$Included, levels = c("Pass", "Fail"))

# rename levels of value
levels(qcMelt$Included) <- c("Included", "Excluded")

with(qcMelt, table(PrimaryDiagnosis, Motion.Exclusion.Level, Included))
## , , Included = Included
## 
##                 Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
##              TD     151     308  372
##              ASD     29     114  173
## 
## , , Included = Excluded
## 
##                 Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
##              TD     221      64    0
##              ASD    144      59    0
qcMelt$Diagnosis4Levels <- factor(qcMelt$Diagnosis4Levels)

Check consistency of diagnosis variables created by Ben and Liwei

#limit to complete predictor cases
dat3 <- filter(dat3, CompletePredictorCases==1)

#convert Diagnosis4Levels to factor
dat3$Diagnosis4Levels <- factor(dat3$Diagnosis4Levels, levels = c("Autism", "Autism+ADHD", "None"))
table(dat3$Diagnosis4Levels)
## 
##      Autism Autism+ADHD        None 
##          49         102         353
table(dat3$ADHD_Secondary)
## 
##   0   1 
## 402 102

49 Autism + 353 None = 402 ADHD_Secondary=0

1. Impact of motion QC on sample size

completeCases <- filter(qcMelt, CompletePredictorCases==1)

tabN<- tableby(PrimaryDiagnosis~Motion.Exclusion.Level, data=filter(completeCases,
                                                                    Included=="Included"))

summary(tabN,  
        title='N',digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)
N
TD (N=791) ASD (N=287)
Motion.Exclusion.Level
   Strict 142 (18.0%) 29 (10.1%)
   Lenient 296 (37.4%) 107 (37.3%)
   None 353 (44.6%) 151 (52.6%)

1b. Summarize sociodemographic info for complete cases for paper

#make M reference level for sex
completeCases$Sex <- relevel(as.factor(completeCases$Sex), "M")

#use chisq for Sex and handedness, kruskal-wallis rank test for Age
tabSex<- tableby( PrimaryDiagnosis~Sex+AgeAtScan+handedness, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))

#use fisher's exact test for Race, k-w for SES
tabRace<- tableby( PrimaryDiagnosis~Race2+SES.Family, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="fe", total=FALSE))

tab12 <- merge(tabSex, tabRace)

#Currently on Stimulants - no test because no TDs are currently on stimulants by design
completeCases$CurrentlyOnStimulants <- factor(completeCases$CurrentlyOnStimulants)

#rename factor levels
levels(completeCases$CurrentlyOnStimulants)[levels(completeCases$CurrentlyOnStimulants)=="0"] <- "No"
levels(completeCases$CurrentlyOnStimulants)[levels(completeCases$CurrentlyOnStimulants)=="1"] <- "Yes"

tabStim<- tableby( PrimaryDiagnosis~CurrentlyOnStimulants, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(total=FALSE, test=FALSE))

tab123 <- merge(tab12, tabStim)
summary(tab123)
TD (N=353) ASD (N=151) p value
Sex < 0.001
   M 245 (69.4%) 127 (84.1%)
   F 108 (30.6%) 24 (15.9%)
AgeAtScan 0.826
   Mean (SD) 10.363 (1.248) 10.324 (1.363)
   Range 8.020 - 12.980 8.010 - 12.990
handedness 0.253
   Left 17 (4.8%) 12 (7.9%)
   Mixed 19 (5.4%) 11 (7.3%)
   Right 317 (89.8%) 128 (84.8%)
Race2 0.004
   African American 36 (10.2%) 9 (6.0%)
   Asian 27 (7.6%) 3 (2.0%)
   Biracial 45 (12.7%) 12 (7.9%)
   Caucasian 245 (69.4%) 127 (84.1%)
SES.Family 0.006
   Mean (SD) 54.135 (9.390) 51.964 (9.379)
   Range 18.500 - 66.000 27.000 - 66.000
CurrentlyOnStimulants
   No 353 (100.0%) 97 (64.2%)
   Yes 0 (0.0%) 54 (35.8%)

Generate version of table to paste into overleaf

tab <- xtable(as.data.frame(summary(tab123)))
print(tab, type="latex")
## % latex table generated in R 4.1.1 by xtable 1.8-4 package
## % Tue Sep 28 21:33:32 2021
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllll}
##   \hline
##  &  & TD (N=353) & ASD (N=151) & p value \\ 
##   \hline
## 1 & **Sex** &  &  & $<$ 0.001 \\ 
##   2 & \&nbsp;\&nbsp;\&nbsp;M & 245 (69.4\%) & 127 (84.1\%) &  \\ 
##   3 & \&nbsp;\&nbsp;\&nbsp;F & 108 (30.6\%) & 24 (15.9\%) &  \\ 
##   4 & **AgeAtScan** &  &  & 0.826 \\ 
##   5 & \&nbsp;\&nbsp;\&nbsp;Mean (SD) & 10.363 (1.248) & 10.324 (1.363) &  \\ 
##   6 & \&nbsp;\&nbsp;\&nbsp;Range & 8.020 - 12.980 & 8.010 - 12.990 &  \\ 
##   7 & **handedness** &  &  & 0.253 \\ 
##   8 & \&nbsp;\&nbsp;\&nbsp;Left & 17 (4.8\%) & 12 (7.9\%) &  \\ 
##   9 & \&nbsp;\&nbsp;\&nbsp;Mixed & 19 (5.4\%) & 11 (7.3\%) &  \\ 
##   10 & \&nbsp;\&nbsp;\&nbsp;Right & 317 (89.8\%) & 128 (84.8\%) &  \\ 
##   11 & **Race2** &  &  & 0.004 \\ 
##   12 & \&nbsp;\&nbsp;\&nbsp;African American & 36 (10.2\%) & 9 (6.0\%) &  \\ 
##   13 & \&nbsp;\&nbsp;\&nbsp;Asian & 27 (7.6\%) & 3 (2.0\%) &  \\ 
##   14 & \&nbsp;\&nbsp;\&nbsp;Biracial & 45 (12.7\%) & 12 (7.9\%) &  \\ 
##   15 & \&nbsp;\&nbsp;\&nbsp;Caucasian & 245 (69.4\%) & 127 (84.1\%) &  \\ 
##   16 & **SES.Family** &  &  & 0.006 \\ 
##   17 & \&nbsp;\&nbsp;\&nbsp;Mean (SD) & 54.135 (9.390) & 51.964 (9.379) &  \\ 
##   18 & \&nbsp;\&nbsp;\&nbsp;Range & 18.500 - 66.000 & 27.000 - 66.000 &  \\ 
##   19 & **CurrentlyOnStimulants** &  &  &  \\ 
##   20 & \&nbsp;\&nbsp;\&nbsp;No & 353 (100.0\%) & 97 (64.2\%) &  \\ 
##   21 & \&nbsp;\&nbsp;\&nbsp;Yes & 0 (0.0\%) & 54 (35.8\%) &  \\ 
##    \hline
## \end{tabular}
## \end{table}

1c. Create exclusion flowchart

1d. Proportion excluded

Define theme for proportion excluded plots

#Easier to make three separate panels and combine them than to use faceting function
My_Theme_prop = theme_light()+theme(
  legend.title =element_blank(),
  axis.title.x = element_text(size = 12),
  axis.title.y = element_text(size = 10),
  plot.title = element_text(size = 30),
  axis.text.x = element_text(size = 8),
  axis.text.y = element_text(size = 8),
  strip.text.x = element_text(size = 10,color="black"),
  strip.background = element_rect(fill="white"))

1d.1 Proportion excluded in each Diagnostic Group

motion <- filter(completeCases, Motion.Exclusion.Level != "None")
motion$Motion.Exclusion.Level <- droplevels(motion$Motion.Exclusion.Level)

levels(motion$Diagnosis4Levels)[levels(motion$Diagnosis4Levels) == "Autism"] = "ASD"
levels(motion$Diagnosis4Levels)[levels(motion$Diagnosis4Levels) == "Autism+ADHD"] = "ASD+ADHD"
levels(motion$Diagnosis4Levels)[levels(motion$Diagnosis4Levels) == "None"] = "TD"

motion$Diagnosis4Levels <- relevel(motion$Diagnosis4Levels, "TD")

motion <- group_by(motion, Diagnosis4Levels, Motion.Exclusion.Level, Included)

dx_proportions <- ggplot(motion, aes(x = Diagnosis4Levels, fill = Included)) + geom_bar(position = "fill",
    alpha = 0.6) + facet_grid(~Motion.Exclusion.Level) + scale_fill_manual(values = c("#FDE599",
    "#9FB0CC")) + scale_color_manual(values = c("#E9D38D", "#8C9AB4")) + My_Theme_prop +
    theme(legend.title = element_blank()) + theme(legend.title = element_blank()) +
    ylab("Proportion of Children") + theme(legend.position = "right") + theme(legend.margin = margin(t = 0,
    r = 0, b = -1, l = -1)) + theme(legend.key.size = unit(0.08, "in")) + theme(axis.title.x = element_blank()) +
    theme(axis.text.x = element_text(size = 10))

png("Figures/fig_propExcludedDx3_cc.png", width = 3, height = 2, units = "in", res = 200)
dx_proportions
invisible(dev.off())

dx_proportions

dx_proportions <- dx_proportions + theme(legend.position = "none") + theme(axis.text.x = element_text(size = 8),
    axis.title.x = element_text(size = 12)) + xlab("Diagnosis Group")

# Pearson's chi squared tests
extib <- tibble(motion)

exNest <- extib %>%
    select(c("Diagnosis4Levels", "Motion.Exclusion.Level", "Included")) %>%
    group_by(Motion.Exclusion.Level) %>%
    tidyr::nest()

# nested models
ex_chisq <- exNest %>%
    mutate(stats = map(data, ~broom::tidy(chisq.test(.x$Diagnosis4Levels, .x$Included)))) %>%
    unnest(stats)

ex_chisq
## # A tibble: 2 × 6
## # Groups:   Motion.Exclusion.Level [2]
##   Motion.Exclusion.Level data               statistic  p.value parameter method 
##   <fct>                  <list>                 <dbl>    <dbl>     <int> <chr>  
## 1 Strict                 <tibble [504 × 2]>      21.2  2.51e-5         2 Pearso…
## 2 Lenient                <tibble [504 × 2]>      12.1  2.34e-3         2 Pearso…

The proportion of children excluded differs across diagnostic groups using both the lenient and strict motion QC.

1d.2 ASD only vs ASD+ADHD

comorbid <- filter(motion, Diagnosis4Levels!="TD")
comorbid$Diagnosis4Levels <- droplevels(comorbid$Diagnosis4Levels)
levels(comorbid$Diagnosis4Levels)
## [1] "ASD"      "ASD+ADHD"
#comorbid <- filter(comorbid, !is.na(Diagnosis4Levels))

coNest <- tibble(comorbid) %>% 
  select(c("Diagnosis4Levels", "CurrentlyOnStimulants",
           "Motion.Exclusion.Level", "Included")) %>% 
  reshape2::melt(id.vars = c("Motion.Exclusion.Level", "Included"),
       measure.vars = c("Diagnosis4Levels", "CurrentlyOnStimulants")) %>% 
  group_by(Motion.Exclusion.Level, variable) %>% 
  tidyr::nest()
## Warning: attributes are not identical across measure variables; they will be
## dropped
#nested models
co_chisq <- coNest %>% 
  mutate(stats = map(data, ~broom::tidy(chisq.test(.x$value, .x$Included)))) %>%
  unnest(stats)

co_chisq
## # A tibble: 4 × 7
## # Groups:   Motion.Exclusion.Level, variable [4]
##   Motion.Exclusion… variable   data    statistic p.value parameter method       
##   <fct>             <fct>      <list>      <dbl>   <dbl>     <int> <chr>        
## 1 Strict            Diagnosis… <tibbl…     0.231  0.631          1 Pearson's Ch…
## 2 Lenient           Diagnosis… <tibbl…     0.463  0.496          1 Pearson's Ch…
## 3 Strict            Currently… <tibbl…     6.40   0.0114         1 Pearson's Ch…
## 4 Lenient           Currently… <tibbl…     6.39   0.0115         1 Pearson's Ch…

1d.3 TD vs ASD only

asdTD <- filter(motion, Diagnosis4Levels!="ASD+ADHD")
asdTD$Diagnosis4Levels <- droplevels(asdTD$Diagnosis4Levels)
levels(asdTD$Diagnosis4Levels)
## [1] "TD"  "ASD"
asdNest <- tibble(asdTD) %>% 
  select(c("Diagnosis4Levels",
           "Motion.Exclusion.Level", "Included")) %>% 
  reshape2::melt(id.vars = c("Motion.Exclusion.Level", "Included"),
       measure.vars = c("Diagnosis4Levels")) %>% 
  group_by(Motion.Exclusion.Level) %>% 
  tidyr::nest()

#nested models
asd_chisq <- asdNest %>% 
  mutate(stats = map(data, ~broom::tidy(chisq.test(.x$value, .x$Included)))) %>%
  unnest(stats)

asd_chisq
## # A tibble: 2 × 6
## # Groups:   Motion.Exclusion.Level [2]
##   Motion.Exclusion.Level data               statistic p.value parameter method  
##   <fct>                  <list>                 <dbl>   <dbl>     <int> <chr>   
## 1 Strict                 <tibble [402 × 3]>      5.04  0.0248         1 Pearson…
## 2 Lenient                <tibble [402 × 3]>      1.56  0.212          1 Pearson…

1d.4 TD vs ASD+ADHD

coTD <- filter(motion, Diagnosis4Levels!="ASD")
coTD$Diagnosis4Levels <- droplevels(coTD$Diagnosis4Levels)
levels(coTD$Diagnosis4Levels)
## [1] "TD"       "ASD+ADHD"
cotNest <- tibble(coTD) %>% 
  select(c("Diagnosis4Levels",
           "Motion.Exclusion.Level", "Included")) %>% 
  reshape2::melt(id.vars = c("Motion.Exclusion.Level", "Included"),
       measure.vars = c("Diagnosis4Levels")) %>% 
  group_by(Motion.Exclusion.Level) %>% 
  tidyr::nest()

#nested models
cot_chisq <- cotNest %>% 
  mutate(stats = map(data, ~broom::tidy(chisq.test(.x$value, .x$Included)))) %>%
  unnest(stats)

cot_chisq
## # A tibble: 2 × 6
## # Groups:   Motion.Exclusion.Level [2]
##   Motion.Exclusion.Level data               statistic  p.value parameter method 
##   <fct>                  <list>                 <dbl>    <dbl>     <int> <chr>  
## 1 Strict                 <tibble [455 × 3]>      16.7  4.33e-5         1 Pearso…
## 2 Lenient                <tibble [455 × 3]>      10.7  1.07e-3         1 Pearso…

1d.5 Proportion excluded by Currently on stimulants in ASD group

stim <- filter(motion, PrimaryDiagnosis=="ASD")
stim$PrimaryDiagnosis <- droplevels(stim$PrimaryDiagnosis)
stim$CurrentlyOnStimulants <- factor(stim$CurrentlyOnStimulants)
levels(stim$CurrentlyOnStimulants)
## [1] "No"  "Yes"
levels(stim$CurrentlyOnStimulants) = c("No", "Yes")

stim <- filter(stim, !is.na(CurrentlyOnStimulants))

stim_proportions <- ggplot(stim, aes(x = CurrentlyOnStimulants, fill = Included))+
  geom_bar(position = "fill", alpha = .6)+
  facet_grid(~Motion.Exclusion.Level)+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme_prop+
  xlab("Currently On Stimulants?")+
  theme(legend.title = element_blank())+
  #theme(legend.position="top")+
  theme(legend.text = element_text(size=7))+
  ylab("Proportion of Autistic Children")+
  theme(legend.position = "none")
  #theme(plot.margin = unit(c(1, 6, 1, 6), "cm"))

prop_legend = cowplot::get_legend(stim_proportions + guides(color = guide_legend(ncol = 1))+
                                    theme(legend.position = "right", 
                                          legend.text = element_text(size = 9),
                                          legend.key.size=unit(.1, "in")))

1d.6 Combine proportion plots into one figure & save

2. rs-fMRI exclusion probability as a function of age and symptom severity

Covariates

phenoVariables <- c("ID", "PrimaryDiagnosis", 
                    "ADOS.Comparable.Total", 
                    "SRS.Score", 
                    "PANESS.TotalOverflowNotAccountingForAge", 
                    "DuPaulHome.InattentionRaw",
                    "DuPaulHome.HyperactivityRaw", 
                    "AgeAtScan", 
                    "WISC.GAI",
                    "Motion.Exclusion.Level", "Included")

phenoIDs <- c("ID", "PrimaryDiagnosis", "Motion.Exclusion.Level", "Included")

aim1 <- reshape2::melt(completeCases[, phenoVariables],
                         id.vars=names(completeCases)[which(names(completeCases) %in% phenoIDs)])

levels(aim1$variable) <- c("ADOS", "SRS", "Motor Overflow", "Inattention", 
                           "Hyperactivity", "Age", "GAI")

aim1G <- group_by(aim1, PrimaryDiagnosis, Motion.Exclusion.Level, Included, variable)

2a. Univarate GAMs

aim1$delta = rep(NA,length=nrow(aim1))
aim1$delta = ifelse(aim1$Included=="Included",1,0)
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level!="None"))
aim1tib$Motion.Exclusion.Level <- droplevels(aim1tib$Motion.Exclusion.Level)

aim1Nest <- aim1tib %>% 
  group_by(Motion.Exclusion.Level, variable) %>% 
  tidyr::nest()

#nested models
nested_gams <- aim1Nest %>% 
  mutate(model = map(data, ~mgcv::gam(1-delta~s(value, k=-1), data = na.omit(.x), 
                                      family=binomial(link=logit), method="REML")), 
         coefs = map(model, tidy, conf.int = FALSE),
         Rsq = map_dbl(model, ~summary(.)$r.sq)) %>% 
  unnest(coefs)

#Ben: correct for 7 lenient and 7 strict 
nested_gams_len <-  nested_gams %>% 
  filter(Motion.Exclusion.Level=="Lenient")

nested_gams_len$p.fdr = p.adjust(nested_gams_len$p.value, method = "BH")

nested_gams_strict <-  nested_gams %>% 
  filter(Motion.Exclusion.Level=="Strict")
  
nested_gams_strict$p.fdr = p.adjust(nested_gams_strict$p.value, method = "BH")

#combine to print
nested_gams <- rbind(nested_gams_len, nested_gams_strict)

#list adjusted p values
nested_gams[, c(1:2,6:11)]
## # A tibble: 14 × 8
## # Groups:   Motion.Exclusion.Level, variable [14]
##    Motion.Exclusion.… variable     edf ref.df statistic  p.value    Rsq    p.fdr
##    <fct>              <fct>      <dbl>  <dbl>     <dbl>    <dbl>  <dbl>    <dbl>
##  1 Lenient            ADOS        2.04   2.47     18.4   2.08e-4 0.0442  4.86e-4
##  2 Lenient            SRS         1.83   2.29     17.4   2.95e-4 0.0466  5.17e-4
##  3 Lenient            Motor Ove…  1.00   1.00     15.6   7.58e-5 0.0308  4.54e-4
##  4 Lenient            Inattenti…  1.38   1.67      7.44  1.17e-2 0.0157  1.17e-2
##  5 Lenient            Hyperacti…  2.15   2.70     13.2   4.59e-3 0.0242  5.36e-3
##  6 Lenient            Age         1.00   1.00      9.33  2.26e-3 0.0164  3.17e-3
##  7 Lenient            GAI         1.00   1.00     14.7   1.30e-4 0.0288  4.54e-4
##  8 Strict             ADOS        1.00   1.00     21.9   2.79e-6 0.0444  9.76e-6
##  9 Strict             SRS         1.00   1.00     22.7   1.53e-6 0.0567  9.76e-6
## 10 Strict             Motor Ove…  1.00   1.00     10.2   1.42e-3 0.0194  1.99e-3
## 11 Strict             Inattenti…  1.00   1.00     17.8   2.47e-5 0.0360  4.31e-5
## 12 Strict             Hyperacti…  1.88   2.35     25.0   1.30e-5 0.0516  3.04e-5
## 13 Strict             Age         1.76   2.20      7.73  2.84e-2 0.0129  2.84e-2
## 14 Strict             GAI         1.00   1.00      7.55  6.02e-3 0.0130  7.02e-3
nested_gams <- nested_gams %>% 
           mutate(LB = map(data, ~round(min(na.omit(.x$value)))),
                  UB = map(data, ~round(max(na.omit(.x$value)))),
                  range = map2(LB, UB, ~seq(from=.x, to=.y, by=1)),
                  logpredict = map2(model, range, ~predict(.x, newdata = data.frame(value = .y), type="link",se.fit=TRUE)),
                  fit = map(logpredict, ~plogis(.x$fit)),
                  lCI = map(logpredict, ~plogis(.x$fit-1.96*.x$se.fit)),
                  hCI = map(logpredict, ~plogis(.x$fit+1.96*.x$se.fit)))

Define theme for gam plots

2a.1 Probability of exclusion as a function of ADOS

ados <- nested_gams %>% 
  filter(variable=="ADOS") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_ados <- ggplot(ados, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='Probability of Exclusion', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("ADOS")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))

2a.2 Probability of exclusion as a function of SRS

srs <- nested_gams %>% 
  filter(variable=="SRS") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_srs <- ggplot(srs, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+
  gam_theme+
  ggtitle("SRS")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.3 Probability of exclusion as a function of Inattention

ina <- nested_gams %>% 
  filter(variable=="Inattention") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_in <- ggplot(ina, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Inattention")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.4 Probability of exclusion as a function of Hyperactivity

hi <- nested_gams %>% 
  filter(variable=="Hyperactivity") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_hi <- ggplot(hi, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Hyperactivity")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.5 Probability of exclusion as a function of Motor Overflow

mo <- nested_gams %>% 
  filter(variable=="Motor Overflow") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_mo <- ggplot(mo, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Motor Overflow")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.6 Probability of exclusion as a function of Age

age<- nested_gams %>% 
  filter(variable=="Age") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_age <- ggplot(age, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Age")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.7 Probability of exclusion as a function of GAI

gai <- nested_gams %>% 
  filter(variable=="GAI") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_gai <- ggplot(gai, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("GAI")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

p_legend = cowplot::get_legend(p_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in")))

2b. Density plots of covariates used to fit GAMs (across included & excluded children)

Define theme for density plots of covariates across included and excluded children

2b.1 ADOS density

ddata <- nested_gams %>% 
  filter(variable=="ADOS") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data) %>% 
  filter(PrimaryDiagnosis=="ASD")

ddata$PrimaryDiagnosis <- droplevels(ddata$PrimaryDiagnosis)

d_ados=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0), limits = c(0, .09), breaks=seq(0, .08, by=.02))+  
  labs(x='', y='Density')+
  scale_fill_manual(values = c("#FDE599"))+ 
  scale_color_manual(values = c("#E9D38D"))+ 
  den_theme

2b.2 SRS density

ddata <- nested_gams %>% 
  filter(variable=="SRS") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_srs=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0), limits=c(0,max(srs$range)),breaks = seq(0, 100 , by = 50))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.3 Inattention density

ddata <- nested_gams %>% 
  filter(variable=="Inattention") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_in=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  labs(x='')+  
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  theme(axis.title.y = element_blank())+
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.4 Hyperactivity/Impulsivity Density

ddata <- nested_gams %>% 
  filter(variable=="Hyperactivity") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_hi=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.5 Motor Overflow Density

ddata <- nested_gams %>% 
  filter(variable=="Motor Overflow") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_mo=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  theme(axis.title.y = element_blank())+
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.6 Age Density

ddata <- nested_gams %>% 
  filter(variable=="Age") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_age=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0), limits=c(8,13), breaks = seq(8, 13 , by = 1))+  
  scale_y_continuous(expand = c(0, 0), limits=c(0,.29), breaks=seq(0, .25, by = .05))+ 
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.7 GAI Density

ddata <- nested_gams %>% 
  filter(variable=="GAI") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_gai=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0), limits = c(0, .035), breaks=seq(0., .03, by=.01))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  labs(x='')+  
  theme(axis.title.y = element_blank())

hist_legend = cowplot::get_legend(d_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in")))

combine gam plots with densities & print

top_row <- cowplot::plot_grid(p_ados, p_srs, p_in, p_hi, p_mo, p_age, p_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
bottom_row <- cowplot::plot_grid(d_ados, d_srs, d_in, d_hi, d_mo, d_age, d_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
## Warning: Removed 91 rows containing non-finite values (stat_density).
## Warning: Removed 19 rows containing non-finite values (stat_density).
png("Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))
dev.off()
## quartz_off_screen 
##                 2
#png("~/Dropbox/Apps/Overleaf/MotionSelectionBias_rsfMRI/Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))

#dev.off()

NOTE: 19 rows = # of participants with missing Motor Overflow scores (imputed for deconfounded group difference)

3. Impact of motion QC on phenotypic representation within diagnostic groups

3a. Split violin plots

My_Theme = theme_light()+theme(
  legend.title = element_blank(),
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  axis.text.x = element_text(size = 5),
  axis.text.y = element_text(size = 8),
  strip.text.x = element_text(size = 12, face = "bold", color="black"),
  strip.text.y = element_text(size = 10, color="black"),
  strip.background = element_rect(fill="white"),
  plot.title = element_text(size = 9, hjust = 0.5))

3a.1 ADOS (ASD only) split violin

NOTE: Missing ADOS scores for one participant evaluated at CARD

ados <- aim1 %>% 
  filter(PrimaryDiagnosis=="ASD" & variable=="ADOS") %>% 
  dplyr::select(-c("PrimaryDiagnosis", "variable"))

ados_violin <- ggplot(ados, aes(Motion.Exclusion.Level, value, fill = Included, color = Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 1.5) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(legend.text = element_text(size = 8))+
  theme(legend.position = "none")+
  ggtitle("ADOS")

3a.2 SRS split violin

srs <- filter(aim1G, variable=="SRS")
aim1_srs <- ggplot(srs, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~., scales = "fixed")+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  # geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("SRS")
  
v_legend = cowplot::get_legend(aim1_srs + guides(color = guide_legend(nrow = 2))+
                                 theme(legend.position = "left", 
                                       legend.text = element_text(size = 10), 
                                       legend.key.size=unit(.1, "in")))
## Warning: Removed 273 rows containing non-finite values (stat_ydensity).
## Warning: Removed 273 rows containing non-finite values (stat_summary).

## Warning: Removed 273 rows containing non-finite values (stat_summary).

NOTE: 273/3 = 91, # of participants missing SRS

with(dat3, table(PrimaryDiagnosis, is.na(SRS.Score)))
##                 
## PrimaryDiagnosis FALSE TRUE
##              TD    263   90
##              ASD   150    1

3a.3 Inattention Split Violin

inat <- filter(aim1G, variable=="Inattention")
paper_inat <- ggplot(inat, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Inattention")

3a.4 Hyperactivity/Impulsivity Spilt Violin

hyp <- filter(aim1G, variable=="Hyperactivity")
paper_hyp <- ggplot(hyp, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Hyperactivity")

3a.5 Motor Overflow Split Violin

overflow <- filter(aim1G, variable=="Motor Overflow")
aim1_of <- ggplot(overflow, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("Motor Overflow")

3a.6 Age Split Violin

age <- filter(aim1G, variable=="Age")
aim1_age <- ggplot(age, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("Age")

3a.7 GAI Split Violin

gai <- filter(aim1G, variable=="GAI")
aim1_gai <- ggplot(gai, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(legend.position = "none")+
  ggtitle("GAI")

Combine split violins into one figure & save

3b. Mann-Whitney U tests to compare included vs excluded participants using lenient motion QC (13 tests)

#run lenient tests first
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Lenient"))

aim1MW <- aim1tib %>% 
  group_by(PrimaryDiagnosis, variable) %>% 
  tidyr::nest()

#hypothesis: included children will have less severe symptoms. NOTE: ADOS only collected in ASD (9 tests)
nested_mw_less <- aim1MW %>% 
  filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
  "Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

#hypothesis: included children will be older and have greater GAI (4 tests)
nested_mw_greater<- aim1MW %>% 
  filter(variable %in% c("Age", "GAI")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

complete_mw <- rbind(nested_mw_less, nested_mw_greater)

complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")

complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:9)]
## # A tibble: 13 × 7
## # Groups:   PrimaryDiagnosis, variable [13]
##    PrimaryDiagnosis variable       statistic  p.value method  alternative  p.fdr
##    <fct>            <fct>              <dbl>    <dbl> <chr>   <chr>        <dbl>
##  1 ASD              ADOS               1780. 0.00921  Wilcox… less        0.0264
##  2 ASD              SRS                1720  0.00579  Wilcox… less        0.0264
##  3 ASD              Motor Overflow     1260  0.000948 Wilcox… less        0.0123
##  4 ASD              Inattention        2442. 0.642    Wilcox… less        0.642 
##  5 ASD              Hyperactivity      2419  0.606    Wilcox… less        0.642 
##  6 ASD              Age                2848  0.0216   Wilcox… greater     0.0352
##  7 ASD              GAI                2960  0.00655  Wilcox… greater     0.0264
##  8 TD               SRS                4562. 0.431    Wilcox… less        0.509 
##  9 TD               Motor Overflow     7086. 0.0569   Wilcox… less        0.0822
## 10 TD               Inattention        8004  0.268    Wilcox… less        0.349 
## 11 TD               Hyperactivity      7018. 0.0198   Wilcox… less        0.0352
## 12 TD               Age               10074. 0.0102   Wilcox… greater     0.0264
## 13 TD               GAI                9882  0.0202   Wilcox… greater     0.0352

3c. Mann-Whitney U tests to compare included vs excluded participants using strict motion QC (13 tests)

#run lenient tests first
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Strict"))

aim1MW <- aim1tib %>% 
  group_by(variable, PrimaryDiagnosis) %>% 
  tidyr::nest()

#hypothesis: included children will have less severe symptoms
nested_mw_less <- aim1MW %>% 
  filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
  "Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

#hypothesis: included children will be older and have greater GAI
nested_mw_greater<- aim1MW %>% 
  filter(variable %in% c("Age", "GAI")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

complete_mw <- rbind(nested_mw_less, nested_mw_greater)

complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")

complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:9)]
## # A tibble: 13 × 7
## # Groups:   PrimaryDiagnosis, variable [13]
##    PrimaryDiagnosis variable       statistic p.value method   alternative  p.fdr
##    <fct>            <fct>              <dbl>   <dbl> <chr>    <chr>        <dbl>
##  1 ASD              ADOS               1350. 0.0236  Wilcoxo… less        0.0767
##  2 ASD              SRS                1312. 0.0176  Wilcoxo… less        0.0763
##  3 ASD              Motor Overflow     1132. 0.0437  Wilcoxo… less        0.0847
##  4 ASD              Inattention        1582  0.189   Wilcoxo… less        0.223 
##  5 ASD              Hyperactivity      1671  0.322   Wilcoxo… less        0.322 
##  6 ASD              Age                2221  0.0165  Wilcoxo… greater     0.0763
##  7 ASD              GAI                1893  0.280   Wilcoxo… greater     0.303 
##  8 TD               SRS                7274. 0.0456  Wilcoxo… less        0.0847
##  9 TD               Motor Overflow    13570. 0.149   Wilcoxo… less        0.194 
## 10 TD               Inattention       14006. 0.147   Wilcoxo… less        0.194 
## 11 TD               Hyperactivity     12199  0.00122 Wilcoxo… less        0.0159
## 12 TD               Age               16656. 0.0374  Wilcoxo… greater     0.0847
## 13 TD               GAI               16378. 0.0685  Wilcoxo… greater     0.111

4. Edgewise functional connectivity as a function of the covariates

Combine partial correlations with melted rs-fMRI usability and covariate info

startEdgeidx = which(names(dat3)=='r.ic1.ic2')
endEdgeidx = which(names(dat3)=='r.ic29.ic30')
names(dat3)[startEdgeidx:endEdgeidx]
##   [1] "r.ic1.ic2"   "r.ic1.ic4"   "r.ic1.ic8"   "r.ic1.ic13"  "r.ic1.ic14" 
##   [6] "r.ic1.ic15"  "r.ic1.ic17"  "r.ic1.ic19"  "r.ic1.ic21"  "r.ic1.ic22" 
##  [11] "r.ic1.ic24"  "r.ic1.ic25"  "r.ic1.ic26"  "r.ic1.ic27"  "r.ic1.ic28" 
##  [16] "r.ic1.ic29"  "r.ic1.ic30"  "r.ic2.ic4"   "r.ic2.ic8"   "r.ic2.ic13" 
##  [21] "r.ic2.ic14"  "r.ic2.ic15"  "r.ic2.ic17"  "r.ic2.ic19"  "r.ic2.ic21" 
##  [26] "r.ic2.ic22"  "r.ic2.ic24"  "r.ic2.ic25"  "r.ic2.ic26"  "r.ic2.ic27" 
##  [31] "r.ic2.ic28"  "r.ic2.ic29"  "r.ic2.ic30"  "r.ic4.ic8"   "r.ic4.ic13" 
##  [36] "r.ic4.ic14"  "r.ic4.ic15"  "r.ic4.ic17"  "r.ic4.ic19"  "r.ic4.ic21" 
##  [41] "r.ic4.ic22"  "r.ic4.ic24"  "r.ic4.ic25"  "r.ic4.ic26"  "r.ic4.ic27" 
##  [46] "r.ic4.ic28"  "r.ic4.ic29"  "r.ic4.ic30"  "r.ic8.ic13"  "r.ic8.ic14" 
##  [51] "r.ic8.ic15"  "r.ic8.ic17"  "r.ic8.ic19"  "r.ic8.ic21"  "r.ic8.ic22" 
##  [56] "r.ic8.ic24"  "r.ic8.ic25"  "r.ic8.ic26"  "r.ic8.ic27"  "r.ic8.ic28" 
##  [61] "r.ic8.ic29"  "r.ic8.ic30"  "r.ic13.ic14" "r.ic13.ic15" "r.ic13.ic17"
##  [66] "r.ic13.ic19" "r.ic13.ic21" "r.ic13.ic22" "r.ic13.ic24" "r.ic13.ic25"
##  [71] "r.ic13.ic26" "r.ic13.ic27" "r.ic13.ic28" "r.ic13.ic29" "r.ic13.ic30"
##  [76] "r.ic14.ic15" "r.ic14.ic17" "r.ic14.ic19" "r.ic14.ic21" "r.ic14.ic22"
##  [81] "r.ic14.ic24" "r.ic14.ic25" "r.ic14.ic26" "r.ic14.ic27" "r.ic14.ic28"
##  [86] "r.ic14.ic29" "r.ic14.ic30" "r.ic15.ic17" "r.ic15.ic19" "r.ic15.ic21"
##  [91] "r.ic15.ic22" "r.ic15.ic24" "r.ic15.ic25" "r.ic15.ic26" "r.ic15.ic27"
##  [96] "r.ic15.ic28" "r.ic15.ic29" "r.ic15.ic30" "r.ic17.ic19" "r.ic17.ic21"
## [101] "r.ic17.ic22" "r.ic17.ic24" "r.ic17.ic25" "r.ic17.ic26" "r.ic17.ic27"
## [106] "r.ic17.ic28" "r.ic17.ic29" "r.ic17.ic30" "r.ic19.ic21" "r.ic19.ic22"
## [111] "r.ic19.ic24" "r.ic19.ic25" "r.ic19.ic26" "r.ic19.ic27" "r.ic19.ic28"
## [116] "r.ic19.ic29" "r.ic19.ic30" "r.ic21.ic22" "r.ic21.ic24" "r.ic21.ic25"
## [121] "r.ic21.ic26" "r.ic21.ic27" "r.ic21.ic28" "r.ic21.ic29" "r.ic21.ic30"
## [126] "r.ic22.ic24" "r.ic22.ic25" "r.ic22.ic26" "r.ic22.ic27" "r.ic22.ic28"
## [131] "r.ic22.ic29" "r.ic22.ic30" "r.ic24.ic25" "r.ic24.ic26" "r.ic24.ic27"
## [136] "r.ic24.ic28" "r.ic24.ic29" "r.ic24.ic30" "r.ic25.ic26" "r.ic25.ic27"
## [141] "r.ic25.ic28" "r.ic25.ic29" "r.ic25.ic30" "r.ic26.ic27" "r.ic26.ic28"
## [146] "r.ic26.ic29" "r.ic26.ic30" "r.ic27.ic28" "r.ic27.ic29" "r.ic27.ic30"
## [151] "r.ic28.ic29" "r.ic28.ic30" "r.ic29.ic30"
signalFC <- dat3[, c(1,startEdgeidx:endEdgeidx)]

dat2 <- merge(aim1, signalFC, all=TRUE, by = "ID")
fcMelt <- reshape2::melt(dat2[,1:160],
                         id.vars=names(dat2)[1:7],
                         variable.name = "edge",
                         value.name = "fc")

4a. Run nested gams

fctib <- tibble(filter(fcMelt, Motion.Exclusion.Level!="None"))
fctib$Motion.Exclusion.Level <- droplevels(fctib$Motion.Exclusion.Level)

fcNest <- fctib %>% 
  filter(Included=="Included") %>% 
  group_by(variable, Motion.Exclusion.Level, edge) %>% 
  tidyr::nest()

#nested models
fc_gams <- fcNest %>% 
  mutate(model = map(data, ~mgcv::gam(fc~s(value), data = na.omit(.x), method="REML")), 
         coefs = map(model, tidy, conf.int = FALSE)) %>% 
  unnest(coefs)

4b. Plot histograms of edgewise p-values from GAMs

Histogram of edgewise p-values for partial correlations as a function of ADOS.

NOTE: TD scores = 0

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="ADOS") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_ados_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  #scale_fill_manual(labels=c('TD','ASD'), values = c("#FDE599", "#9FB0CC"))+
  #scale_color_manual(labels=c('TD','ASD'), values = c("#E9D38D", "#8C9AB4"))+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  ggtitle("ADOS")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  labs(x='p value', y='Count')+
  theme(axis.title.y = element_text(size = 9))+
  theme(axis.title.x = element_text(size = 7))

4b.1 Histogram of edgewise p-values for partial correlations as a function of SRS

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="SRS") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_srs_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("SRS")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.2 Histogram of p-values for partial correlations as a function of inattention

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Inattention") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_in_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Inattention")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.3 Histogram of p-values for partial correlations as a function of hyperactivity

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Hyperactivity") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_hi_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Hyperactivity")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.4 Histogram of p-values for partial correlations as a function of motor overflow

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Motor Overflow") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_mo_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Motor Overflow")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.5 Histogram of p-values for partial correlations as a function of age

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Age") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_age_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Age")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.6 Histogram of p-values for partial correlations as a function of GAI

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="GAI") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_gai_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("GAI")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

hist_p_legend <- cowplot::get_legend(hist_gai_p + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in"), legend.title = element_blank()))

4b.7 Combine histograms and print

fc_hist <- cowplot::plot_grid(hist_ados_p, hist_srs_p, hist_in_p, hist_hi_p, hist_mo_p, hist_age_p, hist_gai_p, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))

#png("~/Dropbox/Apps/Overleaf/MotionSelectionBias_rsfMRI/Figures/fig_hist_fc.png",width=8,height=3,units="in",res=200)
png("Figures/fig_hist_rfc_facet.png",width=8,height=3,units="in",res=200)

cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))
dev.off()
## quartz_off_screen 
##                 2
cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))